# Script conducts statistical analysis and saves results 
# Last Updated: January 15, 2024

# tidyverse (Version 1.3.1)
if((!"tidyverse" %in% installed.packages()) | # If package is not installed 
   packageVersion("tidyverse") != "1.3.1"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/tidyverse/tidyverse_1.3.1.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("tidyverse")
} else { # Otherwise, load package
  library("tidyverse")
}

# PanelMatch (Version 1.0.0)
if((!"PanelMatch" %in% installed.packages()) | # If package is not installed 
   packageVersion("PanelMatch") != "1.0.0"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/PanelMatch/PanelMatch_1.0.0.tar.gz" 
  install.packages(package_url, repos=NULL, type="source")
  library("PanelMatch")
} else { # Otherwise, load package
  library("PanelMatch")
}

# glue (Version 1.6.1)
if((!"glue" %in% installed.packages()) | # If package is not installed 
   packageVersion("glue") != "1.6.1"){ # Or the installed version is not the correct version
  # Install and load correct version of package
  package_url <- "https://cran.r-project.org/src/contrib/Archive/glue/glue_1.6.1.tar.gz"
  install.packages(package_url, repos=NULL, type="source")
  library("glue")
} else { # Otherwise, load package
  library("glue")
}

# Load data
load("data/analysis_data.RData") 

# Set parameters that are constant across all analyses 
pre_treat_periods <- 3 # Number of periods to match covariates on during pre-treatment period
post_treat_periods <- 4 # Number of periods to estimate effect post-withdrawal 
dv <- "reform_n" # Outcome variables
n_iterations <- 5000 # Number of iterations for bootstrapped confidence interval 
set.seed(123) # Set seed

# Define general function to conduct analysis
analysis_function <- function( 
    D, # Name of treatment variable 
    D_placebo_1, # Name of placebo treatment variable led by 1 year 
    D_placebo_2, # Name of placebo treatment variable led by 2 year 
    D_placebo_3, # Name of placebo treatment variable led by 3 year 
    D_placebo_4, # Name of placebo treatment variable led by 4 year 
    D_placebo_5, # Name of placebo treatment variable led by 5 year 
    overlap_var_name, # Name of treatment overlap variable 
    results_file_name # File name of results 
  ){
    # Construct models
    overlap_var_term <- str_c(c("I(lag(", overlap_var_name, ", 1:pre_treat_periods))"), collapse = "")
    balance_var_list <- c("reform_n", "n_members_log", overlap_var_name, 
                          "last_reform", "cumulative_agreements",
                          "pref_mean", "pref_sd",
                          "human_rights", "environment","security", "economics")   
    covariates <- str_c(c("I(lag(n_members_log, 1:pre_treat_periods))", 
                          overlap_var_term,
                          "I(lag(pref_mean, 1:pre_treat_periods))",
                          "I(lag(pref_sd, 1:pre_treat_periods))",
                          "human_rights", "environment", "cumulative_agreements",
                          "security", "economics"), collapse = " + ")
    fmla_true <- as.formula(glue("~ I(lag(reform_n, 1:pre_treat_periods)) + 
                                 {covariates} + I(lag(last_reform, 1:pre_treat_periods))")) 
    
    # Run analyses 
    ## Without refinement
    matches_noref <- PanelMatch(
      exact.match.variables = "excluder",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D, outcome.var = dv,
      refinement.method = "none", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true,
      lead = 0:post_treat_periods, 
      use.diagonal.variance.matrix = T)
    ests_noref <- PanelEstimate(
      sets = matches_noref, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    balance_noref <- get_covariate_balance(
      matches_noref$att, data = analysis_dat, 
      covariates = balance_var_list, plot = F, legend = F) 
    matches_noref_p2 <- PanelMatch(
      exact.match.variables = "excluder_2",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_3, outcome.var = dv, 
      refinement.method = "none", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true, lead = 0:1, 
      use.diagonal.variance.matrix = T)
    ests_noref_p2 <- PanelEstimate(
      sets = matches_noref_p2, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    ## Mahalanobis with up to 5 matches
    matches_maha5 <- PanelMatch(
      exact.match.variables = "excluder",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D, outcome.var = dv,
      refinement.method = "mahalanobis", 
      data = analysis_dat, qoi = "att", size.match = 5,
      covs.formula = fmla_true,
      lead = 0:post_treat_periods, 
      use.diagonal.variance.matrix = T)
    ests_maha5 <- PanelEstimate(
      sets = matches_maha5, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    balance_maha5 <- get_covariate_balance(
      matches_maha5$att, data = analysis_dat,
      covariates = balance_var_list, plot = F, legend = F) 
    matches_maha5_p2 <- PanelMatch(
      exact.match.variables = "excluder_2",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_3, outcome.var = dv,
      refinement.method = "mahalanobis", 
      data = analysis_dat, qoi = "att", size.match = 5,
      covs.formula = fmla_true, lead = 0:1, 
      use.diagonal.variance.matrix = T)
    ests_maha5_p2 <- PanelEstimate(
      sets = matches_maha5_p2, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    ## CBPS Weighting
    matches_cbps <- PanelMatch(
      exact.match.variables = "excluder",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D, outcome.var = dv,
      refinement.method = "CBPS.weight", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true,
      lead = 0:post_treat_periods, 
      use.diagonal.variance.matrix = T)
    ests_cbps <- PanelEstimate(
      sets = matches_cbps, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    balance_cbps <- get_covariate_balance(
      matches_cbps$att, data = analysis_dat,
      covariates = balance_var_list, plot = F, legend = F) 
    matches_cbps_p2 <- PanelMatch(
      exact.match.variables = "excluder_2",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_2, outcome.var = dv,
      refinement.method = "CBPS.weight", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true, lead = 0:1, 
      use.diagonal.variance.matrix = T)
    ests_cbps_p2 <- PanelEstimate(
      sets = matches_cbps_p2, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    ## Supplementary Placebos
    matches_cbps_p1 <- PanelMatch(
      exact.match.variables = "excluder_1",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_1, outcome.var = dv,
      refinement.method = "CBPS.weight", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true, lead = 0, 
      use.diagonal.variance.matrix = T)
    ests_cbps_p1 <- PanelEstimate(
      sets = matches_cbps_p1, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    matches_cbps_p3 <- PanelMatch(
      exact.match.variables = "excluder_3",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_3, outcome.var = dv,
      refinement.method = "CBPS.weight", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true, lead = 0:2, 
      use.diagonal.variance.matrix = T)
    ests_cbps_p3 <- PanelEstimate(
      sets = matches_cbps_p3, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    matches_cbps_p4 <- PanelMatch(
      exact.match.variables = "excluder_4",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_4, outcome.var = dv,
      refinement.method = "CBPS.weight", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true, lead = 0:3, 
      use.diagonal.variance.matrix = T)
    ests_cbps_p4 <- PanelEstimate(
      sets = matches_cbps_p4, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    matches_cbps_p5 <- PanelMatch(
      exact.match.variables = "excluder_5",
      lag = pre_treat_periods, 
      time.id = "year", unit.id = "regime_number_int",
      treatment = D_placebo_5, outcome.var = dv,
      refinement.method = "CBPS.weight", 
      data = analysis_dat, qoi = "att", 
      covs.formula = fmla_true, lead = 0:4, 
      use.diagonal.variance.matrix = T)
    ests_cbps_p5 <- PanelEstimate(
      sets = matches_cbps_p5, data = analysis_dat, 
      number.iterations = n_iterations, confidence.level = .95)
    
    # Save results
    objects_to_save <- c(
      "matches_maha5", "ests_maha5", "balance_maha5", "matches_maha5_p2", "ests_maha5_p2",
      "matches_noref", "ests_noref", "balance_noref", "matches_noref_p2", "ests_noref_p2",
      "matches_cbps", "ests_cbps", "balance_cbps", "matches_cbps_p2", "ests_cbps_p2",
      "matches_cbps_p1", "matches_cbps_p3", "matches_cbps_p4", "matches_cbps_p5",
      "ests_cbps_p1", "ests_cbps_p3", "ests_cbps_p4", "ests_cbps_p5")
    save(list = objects_to_save, file = results_file_name)
}


# Direct Effect  ----

# Define exclusion variable:
# In subsequent analyses, I exclude certain observations from the matched set 
# using exact matching to ensure comparisons between specific groups of regimes. 
# In this first analysis, I allow treated regimes to be matched to any other 
# untreated regime; therefore, the "excluder" is set to a common value across 
# all observations so that all observations matched exactly on this variable.  
analysis_dat$excluder <- 
  analysis_dat$excluder_1 <- 
  analysis_dat$excluder_2 <- 
  analysis_dat$excluder_3 <- 
  analysis_dat$excluder_4 <- 
  analysis_dat$excluder_5 <- 1

# Run analysis 
analysis_function( 
  D = "direct_exit_initial", # Name of treatment variable
  D_placebo_1 = "direct_exit_initial_placebo_1", # Name of placebo treatment variable led by 1 year  
  D_placebo_2 = "direct_exit_initial_placebo_2", # Name of placebo treatment variable led by 2 year  
  D_placebo_3 = "direct_exit_initial_placebo_3", # Name of placebo treatment variable led by 3 year  
  D_placebo_4 = "direct_exit_initial_placebo_4", # Name of placebo treatment variable led by 4 year  
  D_placebo_5 = "direct_exit_initial_placebo_5", # Name of placebo treatment variable led by 5 year 
  overlap_var_name = "overlap", # Name of treatment overlap variable 
  results_file_name = "data/results/results_direct.RData" # File name of results 
)

# Direct Effect (Relative Power == 1) 

# Define exclusion variable:
# Here I define an indicator variable that takes 1 if a regime experiences withdrawal
# by a state that is low relative power and use exact matching to exclude these 
# regimes from the matched set. Same procedure is repeated below. 
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(direct_exit_initial_not_rel_power == 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(direct_exit_initial_placebo_not_rel_power_1 == 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(direct_exit_initial_placebo_not_rel_power_2 == 1, 1, 0),
  excluder_3 = ifelse(direct_exit_initial_placebo_not_rel_power_3 == 1, 1, 0),
  excluder_4 = ifelse(direct_exit_initial_placebo_not_rel_power_4 == 1, 1, 0),
  excluder_5 = ifelse(direct_exit_initial_placebo_not_rel_power_5 == 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_1",  
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_2", 
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_3",  
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_4", 
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_5", 
  overlap_var_name = "overlap_rel_power", 
  results_file_name = "data/results/results_direct_rel_power.RData" 
)

# Direct Effect (Relative Power == 0) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(direct_exit_initial_rel_power == 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(direct_exit_initial_placebo_rel_power_1 == 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(direct_exit_initial_placebo_rel_power_2 == 1, 1, 0),
  excluder_3 = ifelse(direct_exit_initial_placebo_rel_power_3 == 1, 1, 0),
  excluder_4 = ifelse(direct_exit_initial_placebo_rel_power_4 == 1, 1, 0),
  excluder_5 = ifelse(direct_exit_initial_placebo_rel_power_5 == 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_not_rel_power", 
  D_placebo_1 = "direct_exit_initial_placebo_not_rel_power_1",  
  D_placebo_2 = "direct_exit_initial_placebo_not_rel_power_2", 
  D_placebo_3 = "direct_exit_initial_placebo_not_rel_power_3",  
  D_placebo_4 = "direct_exit_initial_placebo_not_rel_power_4", 
  D_placebo_5 = "direct_exit_initial_placebo_not_rel_power_5", 
  overlap_var_name = "overlap_not_rel_power", 
  results_file_name = "data/results/results_direct_not_rel_power.RData" 
)

# Direct Effect (Power Tertile 1) 

# Define exclusion variable
# Here we define an indicator variable that takes 1 if a regime experiences withdrawal
# by a state in another relative power tertile and use exact matching to exclude these 
# regimes from the matched set. Same procedure is repeated below. 
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(direct_exit_initial_rel_power_tertile2 == 1 | direct_exit_initial_rel_power_tertile3 == 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(direct_exit_initial_placebo_rel_power_tertile2_1 == 1 | direct_exit_initial_placebo_rel_power_tertile3_1 == 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(direct_exit_initial_placebo_rel_power_tertile2_2 == 1 | direct_exit_initial_placebo_rel_power_tertile3_2 == 1, 1, 0),
  excluder_3 = ifelse(direct_exit_initial_placebo_rel_power_tertile2_3 == 1 | direct_exit_initial_placebo_rel_power_tertile3_3 == 1, 1, 0),
  excluder_4 = ifelse(direct_exit_initial_placebo_rel_power_tertile2_4 == 1 | direct_exit_initial_placebo_rel_power_tertile3_4 == 1, 1, 0),
  excluder_5 = ifelse(direct_exit_initial_placebo_rel_power_tertile2_5 == 1 | direct_exit_initial_placebo_rel_power_tertile3_5 == 1, 1, 0)
  )

# Run analysis 
analysis_function( 
    D = "direct_exit_initial_rel_power_tertile1", 
    D_placebo_1 = "direct_exit_initial_placebo_rel_power_tertile1_1", 
    D_placebo_2 = "direct_exit_initial_placebo_rel_power_tertile1_2", 
    D_placebo_3 = "direct_exit_initial_placebo_rel_power_tertile1_3", 
    D_placebo_4 = "direct_exit_initial_placebo_rel_power_tertile1_4",
    D_placebo_5 = "direct_exit_initial_placebo_rel_power_tertile1_5", 
    overlap_var_name = "overlap_rel_power_tertile1", 
    results_file_name = "data/results/results_direct_tertile1.RData" 
)

# Direct Effect (Power Tertile 2) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(direct_exit_initial_rel_power_tertile1 == 1 | direct_exit_initial_rel_power_tertile3 == 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_1 == 1 | direct_exit_initial_placebo_rel_power_tertile3_1 == 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_2 == 1 | direct_exit_initial_placebo_rel_power_tertile3_2 == 1, 1, 0),
  excluder_3 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_3 == 1 | direct_exit_initial_placebo_rel_power_tertile3_3 == 1, 1, 0),
  excluder_4 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_4 == 1 | direct_exit_initial_placebo_rel_power_tertile3_4 == 1, 1, 0),
  excluder_5 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_5 == 1 | direct_exit_initial_placebo_rel_power_tertile3_5 == 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power_tertile2", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_tertile2_1",  
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_tertile2_2",  
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_tertile2_3",  
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_tertile2_4",  
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_tertile2_5", 
  overlap_var_name = "overlap_rel_power_tertile2", 
  results_file_name = "data/results/results_direct_tertile2.RData" 
)

# Direct Effect (Power Tertile 3)

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(direct_exit_initial_rel_power_tertile1 == 1 | direct_exit_initial_rel_power_tertile2 == 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_1 == 1 | direct_exit_initial_placebo_rel_power_tertile2_1 == 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_2 == 1 | direct_exit_initial_placebo_rel_power_tertile2_2 == 1, 1, 0),
  excluder_3 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_3 == 1 | direct_exit_initial_placebo_rel_power_tertile2_3 == 1, 1, 0),
  excluder_4 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_4 == 1 | direct_exit_initial_placebo_rel_power_tertile2_4 == 1, 1, 0),
  excluder_5 = ifelse(direct_exit_initial_placebo_rel_power_tertile1_5 == 1 | direct_exit_initial_placebo_rel_power_tertile2_5 == 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power_tertile3", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_tertile3_1",  
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_tertile3_2",  
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_tertile3_3",  
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_tertile3_4",  
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_tertile3_5", 
  overlap_var_name = "overlap_rel_power_tertile3", 
  results_file_name = "data/results/results_direct_tertile3.RData" 
)


# Direct Effect (Restricted) ----

# Define exclusion variable
# Here I exclude indirectly treated regimes from the matched set 
# and replicate the analysis. 
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(indirect_exit_initial == 1 & direct_exit_initial != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(indirect_exit_initial_placebo_1 == 1 & direct_exit_initial_placebo_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(indirect_exit_initial_placebo_2 == 1 & direct_exit_initial_placebo_2 != 1, 1, 0), 
  excluder_3 = ifelse(indirect_exit_initial_placebo_3 == 1 & direct_exit_initial_placebo_3 != 1, 1, 0), 
  excluder_4 = ifelse(indirect_exit_initial_placebo_4 == 1 & direct_exit_initial_placebo_4 != 1, 1, 0), 
  excluder_5 = ifelse(indirect_exit_initial_placebo_5 == 1 & direct_exit_initial_placebo_5 != 1, 1, 0) 
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial", 
  D_placebo_1 = "direct_exit_initial_placebo_1",  
  D_placebo_2 = "direct_exit_initial_placebo_2", 
  D_placebo_3 = "direct_exit_initial_placebo_3", 
  D_placebo_4 = "direct_exit_initial_placebo_4",  
  D_placebo_5 = "direct_exit_initial_placebo_5", 
  overlap_var_name = "overlap", 
  results_file_name = "data/results/results_direct_restricted.RData" 
)

# Direct Effect (Relative Power == 1) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((indirect_exit_initial + direct_exit_initial_not_rel_power > 0) & direct_exit_initial_rel_power != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((indirect_exit_initial_placebo_1 + direct_exit_initial_placebo_not_rel_power_1 > 0) & direct_exit_initial_placebo_rel_power_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((indirect_exit_initial_placebo_2 + direct_exit_initial_placebo_not_rel_power_2 > 0) & direct_exit_initial_placebo_rel_power_2 != 1, 1, 0),
  excluder_3 = ifelse((indirect_exit_initial_placebo_3 + direct_exit_initial_placebo_not_rel_power_3 > 0) & direct_exit_initial_placebo_rel_power_3 != 1, 1, 0),
  excluder_4 = ifelse((indirect_exit_initial_placebo_4 + direct_exit_initial_placebo_not_rel_power_4 > 0) & direct_exit_initial_placebo_rel_power_4 != 1, 1, 0),
  excluder_5 = ifelse((indirect_exit_initial_placebo_5 + direct_exit_initial_placebo_not_rel_power_5 > 0) & direct_exit_initial_placebo_rel_power_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_1",  
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_2", 
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_3",  
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_4", 
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_5", 
  overlap_var_name = "overlap_rel_power", 
  results_file_name = "data/results/results_direct_rel_power_restricted.RData" 
)

# Direct Effect (Relative Power == 0) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((indirect_exit_initial + direct_exit_initial_rel_power > 0) & direct_exit_initial_not_rel_power != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((indirect_exit_initial_placebo_1 + direct_exit_initial_placebo_rel_power_1 > 0) & direct_exit_initial_placebo_not_rel_power_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((indirect_exit_initial_placebo_2 + direct_exit_initial_placebo_rel_power_2 > 0) & direct_exit_initial_placebo_not_rel_power_2 != 1, 1, 0),
  excluder_3 = ifelse((indirect_exit_initial_placebo_3 + direct_exit_initial_placebo_rel_power_3 > 0) & direct_exit_initial_placebo_not_rel_power_3 != 1, 1, 0),
  excluder_4 = ifelse((indirect_exit_initial_placebo_4 + direct_exit_initial_placebo_rel_power_4 > 0) & direct_exit_initial_placebo_not_rel_power_4 != 1, 1, 0),
  excluder_5 = ifelse((indirect_exit_initial_placebo_5 + direct_exit_initial_placebo_rel_power_5 > 0) & direct_exit_initial_placebo_not_rel_power_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_not_rel_power", 
  D_placebo_1 = "direct_exit_initial_placebo_not_rel_power_1",  
  D_placebo_2 = "direct_exit_initial_placebo_not_rel_power_2", 
  D_placebo_3 = "direct_exit_initial_placebo_not_rel_power_3",  
  D_placebo_4 = "direct_exit_initial_placebo_not_rel_power_4", 
  D_placebo_5 = "direct_exit_initial_placebo_not_rel_power_5", 
  overlap_var_name = "overlap_not_rel_power", 
  results_file_name = "data/results/results_direct_not_rel_power_restricted.RData" 
)

# Direct Effect (Power Tertile 1) 

# Define exclusion variable
# Same procedure as before, but now also excluding indirectly treated regimes.
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((indirect_exit_initial == 1 | direct_exit_initial_rel_power_tertile2 == 1 | direct_exit_initial_rel_power_tertile3 == 1) & direct_exit_initial_rel_power_tertile1 != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((indirect_exit_initial_placebo_1 == 1 | direct_exit_initial_placebo_rel_power_tertile2_1 == 1 | direct_exit_initial_placebo_rel_power_tertile3_1 == 1) & direct_exit_initial_placebo_rel_power_tertile1_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((indirect_exit_initial_placebo_2 == 1 | direct_exit_initial_placebo_rel_power_tertile2_2 == 1 | direct_exit_initial_placebo_rel_power_tertile3_2 == 1) & direct_exit_initial_placebo_rel_power_tertile1_2 != 1, 1, 0),
  excluder_3 = ifelse((indirect_exit_initial_placebo_3 == 1 | direct_exit_initial_placebo_rel_power_tertile2_3 == 1 | direct_exit_initial_placebo_rel_power_tertile3_3 == 1) & direct_exit_initial_placebo_rel_power_tertile1_3 != 1, 1, 0),
  excluder_4 = ifelse((indirect_exit_initial_placebo_4 == 1 | direct_exit_initial_placebo_rel_power_tertile2_4 == 1 | direct_exit_initial_placebo_rel_power_tertile3_4 == 1) & direct_exit_initial_placebo_rel_power_tertile1_4 != 1, 1, 0),
  excluder_5 = ifelse((indirect_exit_initial_placebo_5 == 1 | direct_exit_initial_placebo_rel_power_tertile2_5 == 1 | direct_exit_initial_placebo_rel_power_tertile3_5 == 1) & direct_exit_initial_placebo_rel_power_tertile1_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power_tertile1", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_tertile1_1", 
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_tertile1_2", 
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_tertile1_3", 
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_tertile1_4",
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_tertile1_5", 
  overlap_var_name = "overlap_rel_power_tertile1", 
  results_file_name = "data/results/results_direct_tertile1_restricted.RData" 
)

# Direct Effect (Power Tertile 2) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((indirect_exit_initial == 1 | direct_exit_initial_rel_power_tertile1 == 1 | direct_exit_initial_rel_power_tertile3 == 1) & direct_exit_initial_rel_power_tertile2 != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((indirect_exit_initial_placebo_1 == 1 | direct_exit_initial_placebo_rel_power_tertile1_1 == 1 | direct_exit_initial_placebo_rel_power_tertile3_1 == 1) & direct_exit_initial_placebo_rel_power_tertile2_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((indirect_exit_initial_placebo_2 == 1 | direct_exit_initial_placebo_rel_power_tertile1_2 == 1 | direct_exit_initial_placebo_rel_power_tertile3_2 == 1) & direct_exit_initial_placebo_rel_power_tertile2_2 != 1, 1, 0),
  excluder_3 = ifelse((indirect_exit_initial_placebo_3 == 1 | direct_exit_initial_placebo_rel_power_tertile1_3 == 1 | direct_exit_initial_placebo_rel_power_tertile3_3 == 1) & direct_exit_initial_placebo_rel_power_tertile2_3 != 1, 1, 0),
  excluder_4 = ifelse((indirect_exit_initial_placebo_4 == 1 | direct_exit_initial_placebo_rel_power_tertile1_4 == 1 | direct_exit_initial_placebo_rel_power_tertile3_4 == 1) & direct_exit_initial_placebo_rel_power_tertile2_4 != 1, 1, 0),
  excluder_5 = ifelse((indirect_exit_initial_placebo_5 == 1 | direct_exit_initial_placebo_rel_power_tertile1_5 == 1 | direct_exit_initial_placebo_rel_power_tertile3_5 == 1) & direct_exit_initial_placebo_rel_power_tertile2_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power_tertile2", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_tertile2_1",  
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_tertile2_2",  
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_tertile2_3",  
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_tertile2_4",  
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_tertile2_5", 
  overlap_var_name = "overlap_rel_power_tertile2", 
  results_file_name = "data/results/results_direct_tertile2_restricted.RData" 
)

# Direct Effect (Power Tertile 3)

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((indirect_exit_initial == 1 | direct_exit_initial_rel_power_tertile1 == 1 | direct_exit_initial_rel_power_tertile2 == 1) & direct_exit_initial_rel_power_tertile3 != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((indirect_exit_initial_placebo_1 == 1 | direct_exit_initial_placebo_rel_power_tertile1_1 == 1 | direct_exit_initial_placebo_rel_power_tertile2_1 == 1) & direct_exit_initial_placebo_rel_power_tertile3_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((indirect_exit_initial_placebo_2 == 1 | direct_exit_initial_placebo_rel_power_tertile1_2 == 1 | direct_exit_initial_placebo_rel_power_tertile2_2 == 1) & direct_exit_initial_placebo_rel_power_tertile3_2 != 1, 1, 0),
  excluder_3 = ifelse((indirect_exit_initial_placebo_3 == 1 | direct_exit_initial_placebo_rel_power_tertile1_3 == 1 | direct_exit_initial_placebo_rel_power_tertile2_3 == 1) & direct_exit_initial_placebo_rel_power_tertile3_3 != 1, 1, 0),
  excluder_4 = ifelse((indirect_exit_initial_placebo_4 == 1 | direct_exit_initial_placebo_rel_power_tertile1_4 == 1 | direct_exit_initial_placebo_rel_power_tertile2_4 == 1) & direct_exit_initial_placebo_rel_power_tertile3_4 != 1, 1, 0),
  excluder_5 = ifelse((indirect_exit_initial_placebo_5 == 1 | direct_exit_initial_placebo_rel_power_tertile1_5 == 1 | direct_exit_initial_placebo_rel_power_tertile2_5 == 1) & direct_exit_initial_placebo_rel_power_tertile3_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "direct_exit_initial_rel_power_tertile3", 
  D_placebo_1 = "direct_exit_initial_placebo_rel_power_tertile3_1",  
  D_placebo_2 = "direct_exit_initial_placebo_rel_power_tertile3_2",  
  D_placebo_3 = "direct_exit_initial_placebo_rel_power_tertile3_3",  
  D_placebo_4 = "direct_exit_initial_placebo_rel_power_tertile3_4",  
  D_placebo_5 = "direct_exit_initial_placebo_rel_power_tertile3_5", 
  overlap_var_name = "overlap_rel_power_tertile3", 
  results_file_name = "data/results/results_direct_tertile3_restricted.RData" 
)

# Indirect Effect  ----

# Define exclusion variable
# Here I exclude directly treated regimes from the matched set. 
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse(direct_exit_initial == 1 & indirect_exit_initial != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse(direct_exit_initial_placebo_1 == 1 & indirect_exit_initial_placebo_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse(direct_exit_initial_placebo_2 == 1 & indirect_exit_initial_placebo_2 != 1, 1, 0), 
  excluder_3 = ifelse(direct_exit_initial_placebo_3 == 1 & indirect_exit_initial_placebo_3 != 1, 1, 0), 
  excluder_4 = ifelse(direct_exit_initial_placebo_4 == 1 & indirect_exit_initial_placebo_4 != 1, 1, 0), 
  excluder_5 = ifelse(direct_exit_initial_placebo_5 == 1 & indirect_exit_initial_placebo_5 != 1, 1, 0) 
)

# Run analysis 
analysis_function( 
  D = "indirect_exit_initial", 
  D_placebo_1 = "indirect_exit_initial_placebo_1",  
  D_placebo_2 = "indirect_exit_initial_placebo_2", 
  D_placebo_3 = "indirect_exit_initial_placebo_3", 
  D_placebo_4 = "indirect_exit_initial_placebo_4", 
  D_placebo_5 = "indirect_exit_initial_placebo_5", 
  overlap_var_name = "overlap", 
  results_file_name = "data/results/results_indirect.RData" 
)

# Indirect Effect (Absolute Power == 1) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((direct_exit_initial + indirect_exit_initial_not_power > 0) & indirect_exit_initial_power != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((direct_exit_initial_placebo_1 + indirect_exit_initial_placebo_not_power_1 > 0) & indirect_exit_initial_placebo_power_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((direct_exit_initial_placebo_2 + indirect_exit_initial_placebo_not_power_2 > 0) & indirect_exit_initial_placebo_power_2 != 1, 1, 0),
  excluder_3 = ifelse((direct_exit_initial_placebo_3 + indirect_exit_initial_placebo_not_power_3 > 0) & indirect_exit_initial_placebo_power_3 != 1, 1, 0),
  excluder_4 = ifelse((direct_exit_initial_placebo_4 + indirect_exit_initial_placebo_not_power_4 > 0) & indirect_exit_initial_placebo_power_4 != 1, 1, 0),
  excluder_5 = ifelse((direct_exit_initial_placebo_5 + indirect_exit_initial_placebo_not_power_5 > 0) & indirect_exit_initial_placebo_power_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "indirect_exit_initial_power", 
  D_placebo_1 = "indirect_exit_initial_placebo_power_1",  
  D_placebo_2 = "indirect_exit_initial_placebo_power_2", 
  D_placebo_3 = "indirect_exit_initial_placebo_power_3",  
  D_placebo_4 = "indirect_exit_initial_placebo_power_4", 
  D_placebo_5 = "indirect_exit_initial_placebo_power_5", 
  overlap_var_name = "overlap_power", 
  results_file_name = "data/results/results_indirect_power.RData" 
)

# Indirect Effect (Absolute Power == 0) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((direct_exit_initial + indirect_exit_initial_power > 0) & indirect_exit_initial_not_power != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((direct_exit_initial_placebo_1 + indirect_exit_initial_placebo_power_1 > 0) & indirect_exit_initial_placebo_not_power_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((direct_exit_initial_placebo_2 + indirect_exit_initial_placebo_power_2 > 0) & indirect_exit_initial_placebo_not_power_2 != 1, 1, 0),
  excluder_3 = ifelse((direct_exit_initial_placebo_3 + indirect_exit_initial_placebo_power_3 > 0) & indirect_exit_initial_placebo_not_power_3 != 1, 1, 0),
  excluder_4 = ifelse((direct_exit_initial_placebo_4 + indirect_exit_initial_placebo_power_4 > 0) & indirect_exit_initial_placebo_not_power_4 != 1, 1, 0),
  excluder_5 = ifelse((direct_exit_initial_placebo_5 + indirect_exit_initial_placebo_power_5 > 0) & indirect_exit_initial_placebo_not_power_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "indirect_exit_initial_not_power", 
  D_placebo_1 = "indirect_exit_initial_placebo_not_power_1",  
  D_placebo_2 = "indirect_exit_initial_placebo_not_power_2", 
  D_placebo_3 = "indirect_exit_initial_placebo_not_power_3",  
  D_placebo_4 = "indirect_exit_initial_placebo_not_power_4", 
  D_placebo_5 = "indirect_exit_initial_placebo_not_power_5", 
  overlap_var_name = "overlap_not_power", 
  results_file_name = "data/results/results_indirect_not_power.RData" 
)

# Indirect Effect (Power Tertile 1) 

# Define exclusion variable
# Here we define an indicator variable that takes 1 if a regime experiences withdrawal
# by a state in another relative power tertile and use exact matching to exclude these 
# regimes from the matched set. Same procedure is repeated below. 
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((direct_exit_initial + indirect_exit_initial_power_tertile2 + indirect_exit_initial_power_tertile3 > 0) & indirect_exit_initial_power_tertile1 != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile2_1 + indirect_exit_initial_placebo_power_tertile3_1 > 0) & indirect_exit_initial_placebo_power_tertile1_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile2_2 + indirect_exit_initial_placebo_power_tertile3_2 > 0) & indirect_exit_initial_placebo_power_tertile1_2 != 1, 1, 0),
  excluder_3 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile2_3 + indirect_exit_initial_placebo_power_tertile3_3 > 0) & indirect_exit_initial_placebo_power_tertile1_3 != 1, 1, 0),
  excluder_4 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile2_4 + indirect_exit_initial_placebo_power_tertile3_4 > 0) & indirect_exit_initial_placebo_power_tertile1_4 != 1, 1, 0),
  excluder_5 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile2_5 + indirect_exit_initial_placebo_power_tertile3_5 > 0) & indirect_exit_initial_placebo_power_tertile1_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "indirect_exit_initial_power_tertile1", 
  D_placebo_1 = "indirect_exit_initial_placebo_power_tertile1_1", 
  D_placebo_2 = "indirect_exit_initial_placebo_power_tertile1_2", 
  D_placebo_3 = "indirect_exit_initial_placebo_power_tertile1_3", 
  D_placebo_4 = "indirect_exit_initial_placebo_power_tertile1_4",
  D_placebo_5 = "indirect_exit_initial_placebo_power_tertile1_5", 
  overlap_var_name = "overlap_power_tertile1", 
  results_file_name = "data/results/results_indirect_tertile1.RData" 
)

# Indirect Effect (Power Tertile 2) 

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((direct_exit_initial + indirect_exit_initial_power_tertile1 + indirect_exit_initial_power_tertile3 > 0) & indirect_exit_initial_power_tertile2 != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_1 + indirect_exit_initial_placebo_power_tertile3_1 > 0) & indirect_exit_initial_placebo_power_tertile2_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_2 + indirect_exit_initial_placebo_power_tertile3_2 > 0) & indirect_exit_initial_placebo_power_tertile2_2 != 1, 1, 0),
  excluder_3 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_3 + indirect_exit_initial_placebo_power_tertile3_3 > 0) & indirect_exit_initial_placebo_power_tertile2_3 != 1, 1, 0),
  excluder_4 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_4 + indirect_exit_initial_placebo_power_tertile3_4 > 0) & indirect_exit_initial_placebo_power_tertile2_4 != 1, 1, 0),
  excluder_5 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_5 + indirect_exit_initial_placebo_power_tertile3_5 > 0) & indirect_exit_initial_placebo_power_tertile2_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "indirect_exit_initial_power_tertile2", 
  D_placebo_1 = "indirect_exit_initial_placebo_power_tertile2_1",  
  D_placebo_2 = "indirect_exit_initial_placebo_power_tertile2_2",  
  D_placebo_3 = "indirect_exit_initial_placebo_power_tertile2_3",  
  D_placebo_4 = "indirect_exit_initial_placebo_power_tertile2_4",  
  D_placebo_5 = "indirect_exit_initial_placebo_power_tertile2_5", 
  overlap_var_name = "overlap_power_tertile2", 
  results_file_name = "data/results/results_indirect_tertile2.RData" 
)

# Direct Effect (Power Tertile 3)

# Define exclusion variable
analysis_dat <- analysis_dat |> mutate(
  excluder = ifelse((direct_exit_initial + indirect_exit_initial_power_tertile1 + indirect_exit_initial_power_tertile2 > 0) & indirect_exit_initial_power_tertile3 != 1, 1, 0), # Exclude regimes that receive other treatments 
  excluder_1 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_1 + indirect_exit_initial_placebo_power_tertile2_1 > 0) & indirect_exit_initial_placebo_power_tertile3_1 != 1, 1, 0), # Repeat for placebos
  excluder_2 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_2 + indirect_exit_initial_placebo_power_tertile2_2 > 0) & indirect_exit_initial_placebo_power_tertile3_2 != 1, 1, 0),
  excluder_3 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_3 + indirect_exit_initial_placebo_power_tertile2_3 > 0) & indirect_exit_initial_placebo_power_tertile3_3 != 1, 1, 0),
  excluder_4 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_4 + indirect_exit_initial_placebo_power_tertile2_4 > 0) & indirect_exit_initial_placebo_power_tertile3_4 != 1, 1, 0),
  excluder_5 = ifelse((direct_exit_initial + indirect_exit_initial_placebo_power_tertile1_5 + indirect_exit_initial_placebo_power_tertile2_5 > 0) & indirect_exit_initial_placebo_power_tertile3_5 != 1, 1, 0)
)

# Run analysis 
analysis_function( 
  D = "indirect_exit_initial_power_tertile3", 
  D_placebo_1 = "indirect_exit_initial_placebo_power_tertile3_1",  
  D_placebo_2 = "indirect_exit_initial_placebo_power_tertile3_2",  
  D_placebo_3 = "indirect_exit_initial_placebo_power_tertile3_3",  
  D_placebo_4 = "indirect_exit_initial_placebo_power_tertile3_4",  
  D_placebo_5 = "indirect_exit_initial_placebo_power_tertile3_5", 
  overlap_var_name = "overlap_power_tertile3", 
  results_file_name = "data/results/results_indirect_tertile3.RData" 
)

# Load results, give objects specific names, and save 
rm(list = ls()) # Clear environment
load("data/results/results_direct.RData") # Load first batch of results and 
base_var_names <- ls() # Record names of variables
suffix <- "_direct" # Create suffix to identify this set of results
for(i in 1:length(base_var_names)){ # Append suffix to each variable loaded
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_rel_power.RData") # Repeat for other analyses
suffix <- "_direct_rel_power"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_not_rel_power.RData")
suffix <- "_direct_not_rel_power"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_tertile1.RData")
suffix <- "_direct_tertile1"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_tertile2.RData")
suffix <- "_direct_tertile2"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_tertile3.RData")
suffix <- "_direct_tertile3"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_restricted.RData")
suffix <- "_direct_restricted"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_rel_power_restricted.RData") 
suffix <- "_direct_rel_power_restricted"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_not_rel_power_restricted.RData")
suffix <- "_direct_not_rel_power_restricted"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_tertile1_restricted.RData")
suffix <- "_direct_tertile1_restricted"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_tertile2_restricted.RData")
suffix <- "_direct_tertile2_restricted"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_direct_tertile3_restricted.RData")
suffix <- "_direct_tertile3_restricted"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_indirect.RData")
suffix <- "_indirect"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_indirect_power.RData")
suffix <- "_indirect_power"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_indirect_not_power.RData")
suffix <- "_indirect_not_power"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_indirect_tertile1.RData") 
suffix <- "_indirect_tertile1"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_indirect_tertile2.RData") 
suffix <- "_indirect_tertile2"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

load("data/results/results_indirect_tertile3.RData") 
suffix <- "_indirect_tertile3"
for(i in 1:length(base_var_names)){
  var <- eval(parse(text = base_var_names[i]))
  assign(str_c(base_var_names[i], suffix, collapse = ""), var)}

# Save Analyses -----

load("data/analysis_data.RData") # Load data
file_name <- "data/results.RData" # File name of output

objects_to_save <- c("file_name", # List objects to save
                     "analysis_dat", 
                     "unilateral_exits",
                     apropos("ests_"), 
                     apropos("balance_"),
                     apropos("matches_"))
save(list = objects_to_save, file = file_name) # Save
